#include "stdafx.h"

#include "IzhikevichNetwork.h"
#include "PNGData.h"
#include "StreamIO.h"

namespace SpikingNeuronLib
{
	// Public Members -------------------------------------------------------------------------------------------

	array<double>^ IzhikevichNetwork::GetLTPAmplitudes(int neuronIndex)
	{
		array<double>^ result = gcnew array<double>(1001 + _maxDelay);

		// copy LTP values for neuronIndex while adjusting the time by max delay
		for (int i = 0, j = _maxDelay; j < 1001 + _maxDelay; i++, j++)
			result[i] = _LTP[neuronIndex, j];

		return result;
	}

	IzhikevichNetwork::IzhikevichNetwork()
	{
		// defaults:
		//int Ne = 800;			// excitatory neurons
		//int Ni = 200;			// inhibitory neurons
		//int M = 100;			// the number of synapses per neuron
		//int D = 20;			// maximal axonal conduction delay
		unsigned int seed = (unsigned int) DateTime::Now.ToFileTimeUtc();
		this->Init(800, 200, 100, 20, seed);
	}

	IzhikevichNetwork::IzhikevichNetwork(int numExcitatory, int numInhibitory, int numSynapses, int maxDelay)
	{
		unsigned int seed = (unsigned int) DateTime::Now.ToFileTimeUtc();
		this->Init(numExcitatory, numInhibitory, numSynapses, maxDelay, seed);
	}

	IzhikevichNetwork::IzhikevichNetwork(int numExcitatory, int numInhibitory, int numSynapses, int maxDelay, unsigned int seed)
	{
		this->Init(numExcitatory, numInhibitory, numSynapses, maxDelay, seed);
	}

	IzhikevichNetwork::~IzhikevichNetwork()
	{
	}

	void IzhikevichNetwork::SavePostData(String^ filePath)
	{
		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				for (int j = 0; j < _numSynapses; j++)
				{
					sw->Write(_post[i, j].ToString());
					if (j < _numSynapses - 1)
						sw->Write(" ");
					else
						sw->WriteLine();
				}
			}
		}
		finally
		{
			delete sw;
		}
	}

	void IzhikevichNetwork::SaveSynapticData(String^ filePath)
	{
		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				for (int j = 0; j < _numSynapses; j++)
				{
					sw->Write("{0,5:F3}", _s[i, j]);
					if (j < _numSynapses - 1)
						sw->Write(" ");
					else
						sw->WriteLine();
				}
			}
		}
		finally
		{
			delete sw;
		}
	}

	void IzhikevichNetwork::SaveFiringsData(String^ filePath, bool includeOverflow)
	{
		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			for (int i = 0; i < _N_firings; i++)
			{
				// include all lines if includeOverflow is set
				// otherwise include the header [-20 0] and all lines with positive firing time
				if (includeOverflow || _firings[i, 0] == -20 || _firings[i, 0] >= 0)
				{
					sw->WriteLine("{0,3} {1,3}", _firings[i, 0], _firings[i, 1]);
				}
			}
		}
		finally
		{
			delete sw;
		}
	}

	void IzhikevichNetwork::LoadPostData(String^ filePath)
	{
		using namespace System::IO;

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				for (int j = 0; j < _numSynapses; j++)
				{
					Nullable<int> item = StreamIO::ScanStream<int>(sr);
					if (item.HasValue)
					{
						_post[i, j] = item.Value;
					}
					else
					{
						String^ message = String::Format("Error in line length. Expecting {0} items on this line.", _numSynapses);
						throw gcnew ArgumentException(message);
					}
				}
			}
		}
		finally
		{
			delete sr;
		}
	}

	void IzhikevichNetwork::LoadSynapticData(String^ filePath)
	{
		using namespace System::IO;

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				for (int j = 0; j < _numSynapses; j++)
				{
					Nullable<double> item = StreamIO::ScanStream<double>(sr);
					if (item.HasValue)
					{
						_s[i, j] = item.Value;
					}
					else
					{
						String^ message = String::Format("Error in line length. Expecting {0} items on this line.", _numSynapses);
						throw gcnew ArgumentException(message);
					}
				}
			}
		}
		finally
		{
			delete sr;
		}
	}

	void IzhikevichNetwork::LoadFiringsData(String^ filePath)
	{
		using namespace System::IO;

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			int index = 0;
			while ( sr->Peek() >= 0 )
			{
				for (int j = 0; j < 2; j++)
				{
					Nullable<int> item = StreamIO::ScanStream<int>(sr);
					if (item.HasValue)
					{
						_firings[index, j] = item.Value;
					}
					else if (j == 0)
					{
						index--;
						break;	// empty last line
					}
					else
					{
						String^ message = String::Format("Error in line length. Expecting {0} items on this line.", 2);
						throw gcnew ArgumentException(message);
					}
				}
				index++;
			}

			_N_firings = index;	// next available index
		}
		finally
		{
			delete sr;
		}
	}

	void IzhikevichNetwork::SaveNetworkState(String^ filePath)
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// Modifications include porting to a memory-managed environment and refactoring

		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				sw->WriteLine("{0,5:F3} {1,5:F3}", _v[i], _u[i]);
				for (int j = 0; j < _numSynapses; j++)
				{
					sw->WriteLine("{0} {1,5:F3} {2,6:F5}", _post[i, j], _s[i, j], _sd[i, j]);
				}
		
				for (int k = 0; k < _maxDelay; k++)
				{
					sw->Write("{0}  ", _delays_length[i, k]);
					for (int j = 0; j < _delays_length[i, k]; j++)
					{
						sw->Write("{0} ", _delays[i, k, j]);
					}
					sw->WriteLine();
				}

				sw->Write("{0}  ", _N_pre[i]);
				for (int j = 0; j < _N_pre[i]; j++)
				{
					sw->Write("{0} {1} ", _I_pre[i, j], _D_pre[i, j]);
				}
				sw->WriteLine();

				sw->Write(" {0,5:F4} ", _LTD[i]);
				for (int j = 0; j < _maxDelay + 1; j++)
				{
					sw->Write("{0,5:F4} ", _LTP[i, j]);
				}
				sw->WriteLine();
			}

			sw->Write(" {0}", _N_firings);
			for (int i = 0; i < _N_firings; i++)
			{
				sw->Write("{0} {1} ", _firings[i, 0], _firings[i, 1]);
			}
		}
		finally
		{
			delete sw;
		}
	}

	void IzhikevichNetwork::LoadNetworkState(String^ filePath)
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// Modifications include porting to a memory-managed environment, and some bug fixes (see BUGFIX)

		using namespace System::IO;

		// Clear the pre-synaptic arrays
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int j = 0; j < 3 * _numSynapses; j++)
			{
				_I_pre[i, j] = 0;
				_D_pre[i, j] = 0;
				_sd_pre[i, j] = nullptr;
			}
		}

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				_v[i] = StreamIO::ScanStreamForType<double>(sr);
				_u[i] = StreamIO::ScanStreamForType<double>(sr);

				for (int j = 0; j < _numSynapses; j++)
				{
					_post[i, j] = StreamIO::ScanStreamForType<int>(sr);
					_s[i, j] = StreamIO::ScanStreamForType<double>(sr);
					_sd[i, j] = StreamIO::ScanStreamForType<double>(sr);
				}

				for (int k = 0; k < _maxDelay; k++)
				{
					_delays_length[i, k] = StreamIO::ScanStreamForType<int>(sr);

					for (int j = 0; j < _delays_length[i, k]; j++)
					{
						_delays[i, k, j] = StreamIO::ScanStreamForType<int>(sr);
					}
				}

				_N_pre[i] = StreamIO::ScanStreamForType<int>(sr);

				for (int j = 0; j < _N_pre[i]; j++)
				{
					_I_pre[i, j] = StreamIO::ScanStreamForType<int>(sr);
					_D_pre[i, j] = StreamIO::ScanStreamForType<int>(sr);
				}

				_LTD[i] = StreamIO::ScanStreamForType<double>(sr);

				for (int j = 0; j < _maxDelay + 1; j++)
				{
					_LTP[i, j] = StreamIO::ScanStreamForType<double>(sr);
				}
			}
	  
			_N_firings = StreamIO::ScanStreamForType<int>(sr);

			for (int i = 0; i < _N_firings; i++)
			{
				_firings[i, 0] = StreamIO::ScanStreamForType<int>(sr);
				_firings[i, 1] = StreamIO::ScanStreamForType<int>(sr);
			}

			// check for non-truncated firing stack
			// the firing stack is truncated at the end of each network second
			// if a state file has a non-truncated stack, this is evidence of state corruption
			// e.g. a state that has been saved before the clean-up phase
			if (_firings[_N_firings-1, 0] > 0)
			{
				RunWarning(String::Format("Warning: firing stack overflow. The network state in this file may be corrupted: {0}", filePath));
			}
		}
		finally
		{
			delete sr;
		}

		// prepare pre-synaptic references
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int Np = 0; Np < _N_pre[i]; Np++)
			{
				int j = _I_pre[i, Np];
				for (int k = 0; k < _numSynapses; k++)
				{
					if (_post[j, k] == i) 
					{
						// BUGFIX: the Np variable was incorrectly incremented
						//sd_pre[i][Np++]=&sd[j][k];		// Original line
						_sd_pre[i, Np] = gcnew ArrayRef(j, k);
					}
				}
			}
		}
	}

	void IzhikevichNetwork::LoadCrossbarNetworkState(String^ filePath)
	{
		using namespace System::IO;

		// Clear the pre-synaptic arrays
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int j = 0; j < 3 * _numSynapses; j++)
			{
				_I_pre[i, j] = 0;
				_D_pre[i, j] = 0;
				_sd_pre[i, j] = nullptr;
			}
		}

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			Dictionary<int, List<int>^>^ presynapticNeurons = gcnew Dictionary<int, List<int>^>(_totalNeurons);
			Dictionary<int, List<int>^>^ presynapticDelays = gcnew Dictionary<int, List<int>^>(_totalNeurons);
			Dictionary<int, List<int>^>^ delays = gcnew Dictionary<int, List<int>^>(_maxDelay);
			int delaysLength = _numSynapses / _maxDelay;
			int numConnections = StreamIO::ScanStreamForType<int>(sr);

			int offset = 0;
			for (int i = 0; i < _totalNeurons; i++)
			{
				delays->Clear();
				for (int j = 0; j < _numSynapses && offset < numConnections; j++, offset++)
				{
					int preNeuron = StreamIO::ScanStreamForType<int>(sr);
					int postNeuron = StreamIO::ScanStreamForType<int>(sr);
					int delay = StreamIO::ScanStreamForType<int>(sr);
					double weight = StreamIO::ScanStreamForType<double>(sr);
					double derivative = StreamIO::ScanStreamForType<double>(sr);
					int delayIndex = delay - 1;

					System::Diagnostics::Debug::Assert(preNeuron == i);

					_post[preNeuron, j] = postNeuron;
					_s[preNeuron, j] = weight;
					_sd[preNeuron, j] = derivative;

					// track presynaptic neurons (excitatory presynaptic neurons only)
					if (preNeuron < _numExcitatoryNeurons)
					{
						// add to presynaptic neurons
						List<int>^ presynapticNeuronsList;
						if (presynapticNeurons->ContainsKey(postNeuron))
						{
							presynapticNeuronsList = presynapticNeurons[postNeuron];
						}
						else
						{
							presynapticNeuronsList = gcnew List<int>(_numSynapses);
							presynapticNeurons->Add(postNeuron, presynapticNeuronsList);
						}
						presynapticNeuronsList->Add(preNeuron);

						// add to presynaptic delays
						List<int>^ presynapticDelaysList;
						if (presynapticDelays->ContainsKey(postNeuron))
						{
							presynapticDelaysList = presynapticDelays[postNeuron];
						}
						else
						{
							presynapticDelaysList = gcnew List<int>(_numSynapses);
							presynapticDelays->Add(postNeuron, presynapticDelaysList);
						}
						presynapticDelaysList->Add(delayIndex);	// delay (1-20) -> delay index (0-19)
					}

					// track delay offsets for neuron i
					List<int>^ offsets;
					if (delays->ContainsKey(delayIndex))
					{
						offsets = delays[delayIndex];
					}
					else
					{
						offsets = gcnew List<int>(delaysLength);
						delays->Add(delayIndex, offsets);
					}
					offsets->Add(j);
				}

				if (i < _numExcitatoryNeurons)
				{
					// excitatory neurons
					for (int k = 0; k < _maxDelay; k++)
					{
						List<int>^ offsets = delays[k];

						System::Diagnostics::Debug::Assert(offsets->Count == delaysLength);
						_delays_length[i, k] = delaysLength;

						for (int j = 0; j < _delays_length[i, k]; j++)
						{
							_delays[i, k, j] = offsets[j];
						}
					}
				}
				else
				{
					// inhibitory neurons
					List<int>^ offsets = delays[0];

					System::Diagnostics::Debug::Assert(offsets->Count == _numSynapses);
					_delays_length[i, 0] = _numSynapses;

					for (int j = 0; j < _delays_length[i, 0]; j++)
					{
						_delays[i, 0, j] = offsets[j];
					}
				}
			}

			for (int i = 0; i < _totalNeurons; i++)
			{
				_v[i] = StreamIO::ScanStreamForType<double>(sr);
				_u[i] = StreamIO::ScanStreamForType<double>(sr);

				_LTD[i] = StreamIO::ScanStreamForType<double>(sr);

				for (int j = 0; j < _maxDelay + 1; j++)
				{
					_LTP[i, j] = StreamIO::ScanStreamForType<double>(sr);
				}

				List<int>^ presynapticNeuronsList = presynapticNeurons[i];
				List<int>^ presynapticDelaysList = presynapticDelays[i];
				_N_pre[i] = presynapticNeuronsList->Count;

				for (int j = 0; j < _N_pre[i]; j++)
				{
					_I_pre[i, j] = presynapticNeuronsList[j];
					_D_pre[i, j] = presynapticDelaysList[j];
				}
			}
	  
			_N_firings = StreamIO::ScanStreamForType<int>(sr);

			for (int i = 0; i < _N_firings; i++)
			{
				_firings[i, 0] = StreamIO::ScanStreamForType<int>(sr);
				_firings[i, 1] = StreamIO::ScanStreamForType<int>(sr);
				bool isBackgroundFiringEvent = StreamIO::ScanStreamForType<bool>(sr);
			}

			// check for non-truncated firing stack
			// the firing stack is truncated at the end of each network second
			// if a state file has a non-truncated stack, this is evidence of state corruption
			// e.g. a state that has been saved before the clean-up phase
			if (_firings[_N_firings-1, 0] > 0)
			{
				RunWarning(String::Format("Warning: firing stack overflow. The network state in this file may be corrupted: {0}", filePath));
			}
		}
		finally
		{
			delete sr;
		}

		// prepare pre-synaptic references
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int Np = 0; Np < _N_pre[i]; Np++)
			{
				int j = _I_pre[i, Np];
				for (int k = 0; k < _numSynapses; k++)
				{
					if (_post[j, k] == i) 
					{
						// BUGFIX: the Np variable was incorrectly incremented
						//sd_pre[i][Np++]=&sd[j][k];		// Original line
						_sd_pre[i, Np] = gcnew ArrayRef(j, k);
					}
				}
			}
		}
	}

	array<FiringEvent>^ IzhikevichNetwork::GetFiringsData(bool includeOverflow)
	{
		List<FiringEvent>^ events = gcnew List<FiringEvent>();

		// skip over the firing data with times less than 1
		// (unless includeOverflow is set)
		// as this has been copied over from the previous second
		for (int i = 0; i < _N_firings; i++)
		{
			if (includeOverflow || _firings[i, 0] >= 0)
			{
				events->Add(FiringEvent(_firings[i, 0], _firings[i, 1], EventLabel::Unspecified));
			}
		}

		return events->ToArray();
	}

	int IzhikevichNetwork::FindPolychronousGroupsFromTriplet(String^ tag, Triplet^ triplet)
	{
		int maxTime = 0;
		for each (FiringEvent firingEvent in triplet->Events)
		{
			if (firingEvent.Time > maxTime)
				maxTime = firingEvent.Time;
		}

		int presynaptic[W];
		int delays[W];
		for (int i = 0; i < 3; i++)
		{
			FiringEvent firingEvent = triplet->Events[i];
			presynaptic[i] = firingEvent.NeuronIndex;
			// the delays needed in order for the firings of the presynaptic neurons in the triplet
			// to produce convergent input on a postsynaptic target neuron
			delays[i] = maxTime - firingEvent.Time;
		}

		return this->FindPolychronousGroupsFromTriplet(tag, presynaptic, delays);
	}

	int IzhikevichNetwork::FindPolychronousGroupsFromTripletList(String^ tag, List<Triplet^>^ triplets)
	{
		int totalPNGs = 0;

		for each (Triplet^ triplet in triplets)
		{
			totalPNGs += this->FindPolychronousGroupsFromTriplet(tag, triplet);
		}

		return totalPNGs;
	}

	// Polychronous functions--------------------------------
	// The code in this section has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
	// Modifications include:
	//	1. moving related data structures to a memory-managed environment
	//	2. substantial refactoring
	//	3. reporting of PNGs via event mechanism (previously written to console)
	//	4. addition of comments (marked with MG)
	//	5. some bug fixes (see BUGFIX) and additional code that performs bounds checking

	// The original event-triggered (very fast) algorithm
	int IzhikevichNetwork::FindPolychronousGroupsOriginalVersion(String^ tag, int nnum)
	{
		return FindPolychronousGroupsOriginalVersion(tag, C_rel, nnum);
	}

	// The original event-triggered (very fast) algorithm
	int IzhikevichNetwork::FindPolychronousGroupsOriginalVersion(String^ tag, double weightLimit, int nnum)
	{
		int N_polychronous = 0;

		array<int>^ I_my_pre = gcnew array<int>(3 * _numSynapses);
		array<int>^ D_my_pre = gcnew array<int>(3 * _numSynapses);

		// MG get neurons presynaptic to the anchor (nnum) that have 'weighty' synapses
		int N_my_pre = 0;
		for (int i = 0; i < _N_pre[nnum]; i++)
		{
			ArrayRef^ cell = _sd_pre[nnum, i];
			if (_s[cell->X, cell->Y] > weightLimit)
			{
				I_my_pre[N_my_pre] = _I_pre[nnum, i];		// MG presynaptic neuron index
				D_my_pre[N_my_pre] = _D_pre[nnum, i];		// MG delay band
				N_my_pre++;									// MG the number of neurons with weighty synapses
			}
		}

		if (N_my_pre >= W)
		{
			// MG npre is used to create all combinations of triplets
			int npre[W];
			for (int i = 0; i < W; i++)
				npre[i] = i;

			int polylenmax = 1000;				// maximum time range for the group

			array<int>^ group = gcnew array<int>(polylenmax);
			array<int>^ t_fired = gcnew array<int>(polylenmax);
			array<int, 2>^ links = gcnew array<int, 2>(2 * W * polylenmax, 4);

			array<int>^ N_postspikes = gcnew array<int>(polylenmax);
			array<int, 2>^ I_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<int, 2>^ J_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<int, 2>^ D_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<int, 2>^ L_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);

			array<int>^ inh_list = gcnew array<int>(3 * _numInhibitoryNeurons * _numSynapses);
			array<int>^ t_inh_list = gcnew array<int>(3 * _numInhibitoryNeurons * _numSynapses);
			array<int>^ pre_list = gcnew array<int>(3 * _numSynapses);
			array<int>^ del_list = gcnew array<int>(3 * _numSynapses);

   			for (int i = 0; i < polylenmax; i++)
				N_postspikes[i] = 0;

			// MG loop until all combinations of triplets have been explored
			bool searchNextTriplet = true;
			while (searchNextTriplet) 
			{
				this->SearchCurrentTripletOriginalVersion(tag, weightLimit, nnum, npre, N_polychronous, polylenmax,
					I_my_pre, D_my_pre,
					group, t_fired, links,
					N_postspikes, I_postspikes, J_postspikes,
					D_postspikes, L_postspikes,
					inh_list, t_inh_list, pre_list, del_list);

				// MG compute the next combination of available triplets; exit if there are no more
				// MG N_my_pre is the size of the pool of available triplets for the current anchor
				// MG note that combinations can overlap
				int i = 1;
				while (searchNextTriplet && (++npre[W-i] > N_my_pre - i))
				{
					if (++i > W)
					{
						searchNextTriplet = false;		// exit outer while loop
					}
				}

				while (searchNextTriplet && (i > 1))
				{
					npre[W - i + 1] = npre[W - i] + 1;
					i--;
				}
			}

			// set allocated arrays to nullptr
			group = nullptr;
			t_fired = nullptr;
			links = nullptr;
			N_postspikes = nullptr;
			I_postspikes = nullptr;
			J_postspikes = nullptr;
			D_postspikes = nullptr;
			L_postspikes = nullptr;
			inh_list = nullptr;
			t_inh_list = nullptr;
			pre_list = nullptr;
			del_list = nullptr;
		}

		return N_polychronous;
	}

	// The new activity-triggered (slow, but more precise) algorithm
	int IzhikevichNetwork::FindPolychronousGroups(String^ tag, int nnum)
	{
		int N_polychronous = 0;

		array<int>^ I_my_pre = gcnew array<int>(3 * _numSynapses);
		array<int>^ D_my_pre = gcnew array<int>(3 * _numSynapses);

		// MG get neurons presynaptic to the anchor (nnum) that have 'weighty' synapses
		int N_my_pre = 0;
		for (int i = 0; i < _N_pre[nnum]; i++)
		{
			ArrayRef^ cell = _sd_pre[nnum, i];
			if (_s[cell->X, cell->Y] > C_rel)
			{
				I_my_pre[N_my_pre] = _I_pre[nnum, i];		// MG presynaptic neuron index
				D_my_pre[N_my_pre] = _D_pre[nnum, i];		// MG delay band
				N_my_pre++;									// MG the number of neurons with weighty synapses
			}
		}

		// MG there must be at least W weighty synapses for each anchor (W=3)
		// MG candidates are explored in groups of 3 ('triplets')
		if (N_my_pre >= W)
		{
			// MG npre is used to create all combinations of triplets
			int npre[W];
			for (int i = 0; i < W; i++)
				npre[i] = i;

			// MG prepare arrays
			int polylenmax = 1000;				// maximum time range for the group

			array<int>^ group = gcnew array<int>(polylenmax);
			array<int>^ t_fired = gcnew array<int>(polylenmax);
			array<int>^ layer = gcnew array<int>(polylenmax);
			array<int, 2>^ links = gcnew array<int, 2>(2 * W * polylenmax, 4);

			array<int>^ N_postspikes = gcnew array<int>(polylenmax);
			array<int, 2>^ I_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<int, 2>^ J_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<int, 2>^ D_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<int, 2>^ L_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
			array<double, 2>^ C_postspikes = gcnew array<double, 2>(polylenmax, _totalNeurons);

			array<double>^ I = gcnew array<double>(_totalNeurons);

   			for (int i = 0; i < polylenmax; i++)
				N_postspikes[i] = 0;

			// MG loop until all combinations of triplets have been explored
			bool searchNextTriplet = true;
			while (searchNextTriplet) 
			{
				this->SearchCurrentTriplet(tag, nnum, npre, N_polychronous, polylenmax,
					I_my_pre, D_my_pre,
					group, t_fired, layer, links,
					N_postspikes, I_postspikes, J_postspikes,
					D_postspikes, L_postspikes, C_postspikes, I);

				// MG compute the next combination of available triplets; exit if there are no more
				// MG N_my_pre is the size of the pool of available triplets for the current anchor
				// MG note that combinations can overlap
				int i = 1;
				while (searchNextTriplet && (++npre[W-i] > N_my_pre - i))
				{
					if (++i > W)
					{
						searchNextTriplet = false;		// exit outer while loop
					}
				}

				while (searchNextTriplet && (i > 1))
				{
					npre[W - i + 1] = npre[W - i] + 1;
					i--;
				}
			}

			// set allocated arrays to nullptr
			group = nullptr;
			t_fired = nullptr;
			layer = nullptr;
			links = nullptr;
			N_postspikes = nullptr;
			I_postspikes = nullptr;
			J_postspikes = nullptr;
			D_postspikes = nullptr;
			L_postspikes = nullptr;
			C_postspikes = nullptr;
			I = nullptr;
		}

		return N_polychronous;
	}

	int IzhikevichNetwork::FindPolychronousGroupsFromTriplet(String^ tag, int presynaptic[W], int delays[W])
	{
		int N_polychronous = 0;

		int npre[W];
		array<int>^ I_my_pre = gcnew array<int>(3 * _numSynapses);
		array<int>^ D_my_pre = gcnew array<int>(3 * _numSynapses);

		// MG npre is used to index triplet combinations (in this case we have only a single triplet to explore)
		for (int i = 0; i < W; i++)
		{
			npre[i] = i;						// initialise the triplet indexer

			// copy presynaptic neurons and delays
			I_my_pre[i] = presynaptic[i];	// presynaptic neuron index
			D_my_pre[i] = delays[i];		// delay band
		}

		// MG prepare arrays
		int polylenmax = 1000;				// maximum time range for the group

		array<int>^ group = gcnew array<int>(polylenmax);
		array<int>^ t_fired = gcnew array<int>(polylenmax);
		array<int>^ layer = gcnew array<int>(polylenmax);
		array<int, 2>^ links = gcnew array<int, 2>(2 * W * polylenmax, 4);

		array<int>^ N_postspikes = gcnew array<int>(polylenmax);
		array<int, 2>^ I_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
		array<int, 2>^ J_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
		array<int, 2>^ D_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
		array<int, 2>^ L_postspikes = gcnew array<int, 2>(polylenmax, _totalNeurons);
		array<double, 2>^ C_postspikes = gcnew array<double, 2>(polylenmax, _totalNeurons);

		array<double>^ I = gcnew array<double>(_totalNeurons);

   		for (int i = 0; i < polylenmax; i++)
			N_postspikes[i] = 0;

		this->SearchCurrentTriplet(tag, 0, npre, N_polychronous, polylenmax,
			I_my_pre, D_my_pre,
			group, t_fired, layer, links,
			N_postspikes, I_postspikes, J_postspikes,
			D_postspikes, L_postspikes, C_postspikes, I);

		return N_polychronous;
	}

	void IzhikevichNetwork::Shuffle()
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// see BUGFIX for details

		Random^ random = gcnew Random();
		for (int i = 0; i < _numExcitatoryNeurons; i++)
		{
			for (int j = 0; j < _numSynapses; j++)
			{
				if (_post[i, j] < _numExcitatoryNeurons)
				{
					int ri, rj;
					ri = int(System::Math::Floor(random->NextDouble() * _numExcitatoryNeurons));
					do
					{
						rj = int(System::Math::Floor(random->NextDouble() * _numSynapses));
					}
					while (_post[ri, rj] >= _numExcitatoryNeurons);

					// BUGFIX: presumably the intention is to swap random locations; this is not doing as intended
					// original code commented out...
					//double x = _s[ri, rj];
					//double y = _sd[ri, rj];

					//_s[i, j] = _s[ri, rj];
					//_sd[i, j] = _sd[ri, rj];

					//_s[ri, rj] = x;
					//_sd[ri, rj] = y;

					// ...rewritten as a swap
					double x = _s[ri, rj];
					double y = _sd[ri, rj];

					_s[ri, rj] = _s[i, j];
					_sd[ri, rj] = _sd[i, j];

					_s[i, j] = x;
					_sd[i, j] = y;
				}
			}
		}
	}

	void IzhikevichNetwork::RandomizeWeights()
	{
		unsigned int seed = (unsigned int) DateTime::Now.ToFileTimeUtc();
		Random^ random = gcnew Random(seed);

		for (int i = 0; i < _numExcitatoryNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
				_s[i, j] = 10 * random->NextDouble();

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
				_s[i, j] = -(10 * random->NextDouble());

		random = nullptr;
	}

	// Randomize synaptic weights by perturbing the base weight up or down by random values within the specified range
	void IzhikevichNetwork::RandomizeWeights(int range)
	{
		unsigned int seed = (unsigned int) DateTime::Now.ToFileTimeUtc();
		Random^ random = gcnew Random(seed);

		for (int i = 0; i < _numExcitatoryNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
			{
				double direction = random->NextDouble() > 0.5? -1.0 : 1.0;
				double perturbation = range * random->NextDouble();
				_s[i, j] = 6 + (direction * perturbation);
			}

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
			{
				double direction = random->NextDouble() > 0.5? -1.0 : 1.0;
				double perturbation = range * random->NextDouble();
				_s[i, j] = -5 + (direction * perturbation);
			}

		random = nullptr;
	}

	List<Connection^>^ IzhikevichNetwork::GetConnections()
	{
		List<Connection^>^ connections = gcnew List<Connection^>(_totalNeurons * _numSynapses);

		for (int preNeuron = 0; preNeuron < _totalNeurons; preNeuron++)
		{
			for (int delay = 0; delay < _maxDelay; delay++)
			{
				for (int j = 0; j < _delays_length[preNeuron, delay]; j++)
				{
					int offset = _delays[preNeuron, delay, j];
					int postNeuron = _post[preNeuron, offset];
					double weight = _s[preNeuron, offset];

					// map delays 0-19 to delays 1-20
					Connection^ connection = gcnew Connection(preNeuron, postNeuron, delay + 1, weight);
					connections->Add(connection);
				}
			}
		}

		return connections;
	}

	// Private Members -------------------------------------------------------------------------------------------

	void IzhikevichNetwork::Init(int numExcitatory, int numInhibitory, int numSynapses, int maxDelay, unsigned int seed)
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// Modifications include porting to a memory-managed environment and refactoring

		_initialized = false;
		_N_firings_max = 100000;

		Random^ random = gcnew Random(seed);

		_numExcitatoryNeurons = numExcitatory;								// excitatory neurons
		_numInhibitoryNeurons = numInhibitory;								// inhibitory neurons
		_totalNeurons = _numExcitatoryNeurons + _numInhibitoryNeurons;		// total number
		_numSynapses = numSynapses;											// the number of synapses per neuron
		_maxDelay = maxDelay;												// maximal axonal conduction delay

		_a = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _numExcitatoryNeurons; i++)
			_a[i] = 0.02;

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			_a[i] = 0.1;

		_d = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _numExcitatoryNeurons; i++)
			_d[i] = 8;

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			_d[i] = 2;

		// post-synaptic neuron indices
		_post = gcnew array<int, 2>(_totalNeurons, _numSynapses);
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int j = 0; j < _numSynapses; j++) 
			{
				int r;
				bool exists;
				do
				{
					exists = false;
					if (i < _numExcitatoryNeurons)
					{
						// excitatory neurons can connect to any other neuron
						r = random->Next(_totalNeurons);		// 0 to _totalNeurons-1
					}
					else
					{
						// inhibitory neurons can only connect to excitatory neurons
						r = random->Next(_numExcitatoryNeurons);		// 0 to _numExcitatoryNeurons-1
					}

					if (r == i)
					{
						exists = true;
						continue;
					}

					for (int k = 0; k < j; k++)
					{
						if (_post[i, k] == r)
						{
							exists = true;
							break;
						}
					}
				}
				while (exists);

				_post[i, j] = r;
			}
		}

		// initial synaptic weights
		_s = gcnew array<double, 2>(_totalNeurons, _numSynapses);
		for (int i = 0; i < _numExcitatoryNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
				_s[i, j] = 6; 

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
				_s[i, j] = -5;

		// derivatives of synaptic weights
		_sd = gcnew array<double, 2>(_totalNeurons, _numSynapses);
  		for (int i = 0; i < _totalNeurons; i++)
			for (int j = 0; j < _numSynapses; j++)
				_sd[i, j] = 0;

		// axonal delays
		// delays[i][j][k]: the kth offset for neuron i having delay j (where offset is the current synaptic offset into post, s, sd)
		_delays_length = gcnew array<short, 2>(_totalNeurons, _maxDelay);
		_delays = gcnew array<short, 3>(_totalNeurons, _maxDelay, _numSynapses);	// Note: longer in the 3rd dimension than it needs to be
  		for (int i = 0; i < _totalNeurons; i++)
		{
			short ind = 0;
			if (i < _numExcitatoryNeurons)
			{
				for (int j = 0; j < _maxDelay; j++) 
				{
					_delays_length[i, j] = _numSynapses / _maxDelay;
					for (int k = 0; k < _delays_length[i, j]; k++)
						_delays[i, j, k] = ind++;
				}
			}
			else
			{
				for (int j = 0; j < _maxDelay; j++)
					_delays_length[i, j] = 0;

				_delays_length[i, 0] = _numSynapses;

				for (int k = 0; k < _delays_length[i, 0]; k++)
						_delays[i, 0, k] = ind++;
			}
		}
	
		// indices of pre-synaptic excitatory neurons
		// Note: the allocation size of these arrays assumes _N_pre[i] < 3 * _numSynapses
		// LoadNetworkState() also makes this assumption (when clearing these arrays)
		_N_pre = gcnew array<int>(_totalNeurons);
		_I_pre = gcnew array<int, 2>(_totalNeurons, 3 * _numSynapses);
		_D_pre = gcnew array<int, 2>(_totalNeurons, 3 * _numSynapses);
		_sd_pre = gcnew array<ArrayRef^, 2>(_totalNeurons, 3 * _numSynapses);
  		for (int i = 0; i < _totalNeurons; i++)
		{
			_N_pre[i] = 0;

			for (int j = 0; j < _numExcitatoryNeurons; j++)
			{
				for (int k = 0; k < _numSynapses; k++)	// current synaptic offset
				{
					if (_post[j, k] == i) 
					{
						_I_pre[i, _N_pre[i]] = j;	// the indices of neurons that are pre-synaptic to i

						// find the delay, dd, for which _delays[j, dd, jj] produces the current synaptic offset
						for (int dd = 0; dd < _maxDelay; dd++)
						{
							for (int jj = 0; jj < _delays_length[j, dd]; jj++)
							{
								if (_post[j, _delays[j, dd, jj]] == i)
									_D_pre[i, _N_pre[i]] = dd;
							}
						}

						_sd_pre[i, _N_pre[i]++] = gcnew ArrayRef(j, k);	// &_sd[j, k];
					}
				}
			}
		}

		// LTP
		_LTP = gcnew array<double, 2>(_totalNeurons, 1001 + _maxDelay);
		for (int i = 0; i < _totalNeurons; i++)
			for (int j = 0; j < 1 + _maxDelay; j++)
				_LTP[i, j] = 0;

		// LTD
		_LTD = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _totalNeurons; i++)
			_LTD[i] = 0;

		// initial values for v
		_v = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _totalNeurons; i++)
			_v[i] = -65 + 10 * random->NextDouble();		// -65.0 inclusive to -55.0 exclusive

		// initial values for u
		_u = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _totalNeurons; i++)
			_u[i] = 0.2 * _v[i];

		// spike timings
		_firings = gcnew array<int, 2>(_N_firings_max, 2);
		_N_firings = 1;
		_firings[0, 0] = -(_maxDelay);
		_firings[0, 1] = 0;

		random = nullptr;

		_initialized = true;
	}

	int IzhikevichNetwork::RunImpl(int totalSeconds, PatternGenerator^ backgroundPattern, IStimulus^ inputPattern)
	{
		return RunNetwork(0, totalSeconds, backgroundPattern, inputPattern);
	}

	int IzhikevichNetwork::RunImpl(int startSeconds, int endSeconds, PatternGenerator^ backgroundPattern, IStimulus^ inputPattern)
	{
		return RunNetwork(startSeconds, endSeconds, backgroundPattern, inputPattern);
	}

	int IzhikevichNetwork::RunNetwork(int startSeconds, int endSeconds, PatternGenerator^ backgroundPattern, IStimulus^ inputPattern)
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// Modifications include porting to a memory-managed environment and refactoring

		int result = 1;

		if (_initialized)
		{
			array<double>^ I = gcnew array<double>(_totalNeurons);

			for (int sec = startSeconds; sec < endSeconds; sec++)
			{
				// get the external stimulus firing times for the current second (both input and background)
				_stimulusFiringData = nullptr;
				if (inputPattern)
				{
					_stimulusFiringData = inputPattern->GetNeuronFiringData(sec);
				}
				_backgroundFiringData = nullptr;
				if (backgroundPattern)
				{
					_backgroundFiringData = backgroundPattern->GetNeuronFiringData(sec);
				}

				for (int t = 0; t < 1000; t++)		// simulation of 1 sec
				{
					// input vector initialisation
					for (int i = 0; i < _totalNeurons; i++)
						I[i] = 0;

					for (int i = 0; i < _totalNeurons; i++)
					{
						if (_v[i] >= 30)
						{
							_v[i] = -65;
							_u[i] += _d[i];

							_LTP[i, t + _maxDelay] = 0.1;
							_LTD[i] = 0.12;

							for (int j = 0; j < _N_pre[i]; j++)
							{
								//*sd_pre[i][j] += LTP[I_pre[i][j]][t+D-D_pre[i][j]-1];
								ArrayRef^ cell = _sd_pre[i, j];
								_sd[cell->X, cell->Y] += _LTP[_I_pre[i, j], t + _maxDelay - _D_pre[i, j] - 1];
							}

							_firings[_N_firings, 0] = t;
							_firings[_N_firings++, 1] = i;
							if (_N_firings == _N_firings_max)
							{
								RunWarning(String::Format("Two many spikes. Truncating at t={0}", t));
								_N_firings = 1;
							}
						}
					}

					int k = _N_firings - 1;

					while (t - _firings[k, 0] < _maxDelay)
					{
						int currentTime = _firings[k, 0];
						int currentNeuron = _firings[k, 1];
						for (int j = 0; j < _delays_length[currentNeuron, t - currentTime]; j++)
						{
							int currentOffset = _delays[currentNeuron, t - currentTime, j];	// the jth offset for currentNeuron having delay t - currentTime

							int i = _post[currentNeuron, currentOffset];

							I[i] += _s[currentNeuron, currentOffset];		// internal input

							if (currentNeuron < _numExcitatoryNeurons)
							{
								_sd[currentNeuron, currentOffset] -= _LTD[i];
							}
						}

						k--;
					}

					// set external input (internal input has already been set)
					if (_stimulusFiringData)
					{
						// set input vector values (corresponding to specific neurons) for the current millisecond
						if (_stimulusFiringData[t])
						{
							array<int>^ firedNeuron = _stimulusFiringData[t];
							for each (int neuron in firedNeuron)
							{
								// reset _u[neuron] to default so that the neuron will fire even if it has recently fired
								_u[neuron] = -13.0;

								// add whatever value is required to get _v[neuron] to +20
								double current = (0.04 * _v[neuron] + 5) * _v[neuron] + 140 - _u[neuron] + I[neuron];
								if (current < 20)
									I[neuron] += (20 - current);		// MGSC
							}
						}
					}

					// optionally include random background firing
					if (_backgroundFiringData)
					{
						if (_backgroundFiringData[t])
						{
							array<int>^ selectedNeurons = _backgroundFiringData[t];
							for each (int neuron in selectedNeurons)
							{
								I[neuron] += 20;
							}
						}
					}

					for (int i = 0; i < _totalNeurons; i++)
					{
						_v[i] += 0.5 * ((0.04 * _v[i] + 5) * _v[i] + 140 - _u[i] + I[i]);	// for numerical stability
						_v[i] += 0.5 * ((0.04 * _v[i] + 5) * _v[i] + 140 - _u[i] + I[i]);	// time step is 0.5 ms

						_u[i] += _a[i] * (0.2 * _v[i] - _u[i]);

						_LTP[i, t + _maxDelay + 1] = 0.95 * _LTP[i, t + _maxDelay];			// tau = 20 ms
						_LTD[i] *= 0.95;													// tau = 20 ms
					}
				}	// end of loop over 1000 msecs

				// average firing rate for excitatory neurons
				int excitatoryFiringCount = 0;
				for (int i = 1; i < _N_firings; i++)
				{
					if ((_firings[i, 0] >= 0) && (_firings[i, 1] < _numExcitatoryNeurons))
						excitatoryFiringCount++;
				}
				double frate = (double)excitatoryFiringCount / _numExcitatoryNeurons;

				// // proportion of excitatory-excitatory connections that are near maximal synaptic weight
				int tstrength = 0;
				for (int i = 0; i < _numExcitatoryNeurons; i++)
				{
					for (int j = 0; j < _numSynapses; j++)
						if ((_s[i, j] > 0.9 * C_max) && (_post[i, j] < _numExcitatoryNeurons))
							tstrength++;
				}
				double strength = 100 * (double)tstrength / _numExcitatoryNeurons / _numSynapses;

				// average firing rate for inhibitory neurons
				int inhibitoryFiringCount = 0;
				for (int i = 1; i < _N_firings; i++)
				{
					if ((_firings[i, 0] >= 0) && (_firings[i, 1] >= _numExcitatoryNeurons))
						inhibitoryFiringCount++;
				}
				double ifrate = (double)inhibitoryFiringCount / _numInhibitoryNeurons;

				OneSecondTickBeforeCleanup(sec, frate, strength, ifrate);

				if (backgroundPattern)
					backgroundPattern->Reset();		// prepare background firing patterns for the next second

				for (int i = 0; i < _totalNeurons; i++)
					for (int j = 0; j < _maxDelay + 1; j++)
						_LTP[i, j] = _LTP[i, 1000 + j];

				int k = _N_firings - 1;
				while (1000 - _firings[k, 0] < _maxDelay)
				{
					k--;
				}

				for (int i = 1; i < _N_firings - k; i++)
				{
					_firings[i, 0] = _firings[k + i, 0] - 1000;
					_firings[i, 1] = _firings[k + i, 1];
				}
				_N_firings = _N_firings - k;

				for (int i = 0; i < _numExcitatoryNeurons; i++)
					for (int j = 0; j < _numSynapses; j++)
					{
						_sd[i, j] *= 0.9;				// tau = 250 ms
						_s[i, j] += 0.01 + _sd[i, j];
						if (_s[i, j] > C_max)
							_s[i, j] = C_max;
						if (_s[i, j] < 0)
							_s[i, j] = 0;
					}

				// SaveNetworkState() and SaveSynapticData() can be called from within this callback, or at any other time
				// (except from within the beforeCleanupStatusCallback as the above cleanup code must be called first)
				OneSecondTickAfterCleanup(sec);

			}	// end loop over many seconds

			I = nullptr;

			result = 0;
		}

		return result;
	}

	void IzhikevichNetwork::SearchCurrentTripletOriginalVersion(String^ tag, double weightLimit, int nnum, int npre[W], int &N_polychronous, int polylenmax,
			array<int>^ I_my_pre, array<int>^ D_my_pre,
			array<int>^ group, array<int>^ t_fired, array<int, 2>^ links,
			array<int>^ N_postspikes, array<int, 2>^ I_postspikes, array<int, 2>^ J_postspikes,
			array<int, 2>^ D_postspikes, array<int, 2>^ L_postspikes,
			array<int>^ inh_list, array<int>^ t_inh_list, array<int>^ pre_list, array<int>^ del_list)
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// see BUGFIX for details

		const int V = 2;				// number of spikes needed to fire a cell
		const int latency = 4;			// latency from the moment spikes arrive to the moment the neuron fires
		const int tf_final = 2;
		const int inh_span = 5;
		const int slacke = 0;
		const int slacki = 0;

		int timing;

		int Dmax = 0;
		for (int i = 0; i < W; i++)
		{
			if (Dmax < D_my_pre[npre[i]])
				Dmax = D_my_pre[npre[i]];
		}

		for (int i = 0; i < W; i++)
		{
			group[i] = I_my_pre[npre[i]];
			t_fired[i] = Dmax - D_my_pre[npre[i]];
			
			for (int d = 0; d < _maxDelay; d++)
			{
				for (int j = 0; j < _delays_length[group[i], d]; j++)
				{
					int offset = _delays[group[i], d, j];
					int p = _post[group[i], offset];
					if (((_s[group[i], offset] > weightLimit) && (d >= D_my_pre[npre[i]])) || (p >= _numExcitatoryNeurons))
					{
						timing = t_fired[i] + d + 1;
						J_postspikes[timing, N_postspikes[timing]] = group[i];				// presynaptic
						D_postspikes[timing, N_postspikes[timing]] = d;						// delay
						L_postspikes[timing, N_postspikes[timing]] = 1;						// layer
						I_postspikes[timing, N_postspikes[timing]++] = p;					// index of post target	
					}
				}
			}
		}

		int N_inh_list = 0;
		int NL_max = 1;
		int N_links = 0;
		int N_fired = W;
		int t_last = _maxDelay + _maxDelay + latency + 1;

		for (int t = 0; t < t_last; t++)
		{
			while (N_postspikes[t] > 0)
			{
				N_postspikes[t]--;

				int already = 0;

				if (I_postspikes[t, N_postspikes[t]] < _numExcitatoryNeurons)
				{
					for (int j = 0; j < N_fired; j++)
					{
						if ((I_postspikes[t, N_postspikes[t]] == group[j]) && (t_fired[j] > t - 20))
							already = 2;
					}

					for (int j = 0; j < N_inh_list; j++)
					{
						if ((inh_list[j] == I_postspikes[t, N_postspikes[t]]) && (t_inh_list[j] <= t) && (t_inh_list[j] + inh_span >= t))
							already++;
					}
				}
				else
				{
					for (int j = 0; j < N_fired; j++)
					{
						if ((I_postspikes[t, N_postspikes[t]] == group[j]) && (t_fired[j] > t - 20))
							already = 2;
					}
				}
 
				if ((already <= 0) && (N_fired < polylenmax))
				{
					int last;
					int NL = 0;
					int exc = 0;
					int first = -1;
					int t_fire = t + 4;

					for (int tf = 0; tf <= tf_final; tf++)
					{
						for (int j = 0; j <= N_postspikes[t+tf] - 0.5 * tf; j++)
						{
							if (I_postspikes[t + tf, j] == I_postspikes[t, N_postspikes[t]])
							{
								if (first < 0)
									first = tf;
								last = tf;
								pre_list[exc] = J_postspikes[t+tf, j];
								del_list[exc] = D_postspikes[t+tf, j];
								exc++;
								t_fire = t + tf;
								if (NL < L_postspikes[t+tf, j])
									NL = L_postspikes[t+tf, j];
							}
						}
					}

					if (((first + 1 >= last) && (exc >= V - slacke))
						|| (((I_postspikes[t, N_postspikes[t]] < _numExcitatoryNeurons) && (exc >= V))
						|| ((I_postspikes[t, N_postspikes[t]] >= _numExcitatoryNeurons) && (first + 2 >= last) && (exc >= V - slacki))))
					{
						int i = N_fired++;
	 					group[i] = I_postspikes[t, N_postspikes[t]];
						t_fired[i] = t_fire + latency;
						if (exc == 2)
							t_fired[i] = t_fired[i] + latency;		// longer latencies for weaker inputs

						for (int j = 0; j < exc; j++)
						{
							links[N_links, 0] = pre_list[j];
							links[N_links, 1] = group[i];
							links[N_links, 2] = del_list[j];
							links[N_links++, 3] = NL;
						}

						if (group[i] < _numExcitatoryNeurons)
						{
							for (int d = 0; d < _maxDelay; d++)
							{
								for (int j = 0; j < _delays_length[group[i], d]; j++)
								{
									int offset = _delays[group[i], d, j];
									if (_s[group[i], offset] > weightLimit)
									{
										timing = t_fired[i] + d + 1;
										// BUGFIX: it is possible for 'timing' to exceed the maximum index of the _postspikes arrays
										// the next line checks this limit
										if ((timing < polylenmax) && (N_postspikes[timing] < _totalNeurons))	// check the limit
										{
											J_postspikes[timing, N_postspikes[timing]] = group[i];					// presynaptic
											D_postspikes[timing, N_postspikes[timing]] = d;							// delay
											L_postspikes[timing, N_postspikes[timing]] = NL + 1;					// layer
											I_postspikes[timing, N_postspikes[timing]++] = _post[group[i], offset];	// index of post target
										}
									}
								}
							}
						}
						else
						{
							for (int d = 0; d < _maxDelay; d++)
							{
								for (int j = 0; j < _delays_length[group[i], d]; j++)
								{
									inh_list[N_inh_list] = _post[group[i], _delays[group[i], d, j]];
									t_inh_list[N_inh_list++] = t_fired[i] + d;
								}
							}

							int j = 0;
							while (t_inh_list[j] + inh_span < t)
								j++;

							if (j > 10)	// needs cleaning
							{
								for (int k = j; k < N_inh_list; k++)
								{
									inh_list[k - j] = inh_list[k];
									t_inh_list[k - j] = t_inh_list[k];
								}
								N_inh_list = N_inh_list - j;
							}
						}
					
						if (NL > NL_max)
							NL_max = NL;
				
						if (t_last < timing + 1) 
						{
							t_last = timing + 1;
							if (t_last > polylenmax - _maxDelay - latency - tf_final - 1)
								t_last = polylenmax - _maxDelay - latency - tf_final - 1;
						}
					}
				}
			}
		}

		if (t_last == polylenmax - _maxDelay)
		{
			for (int d = t_last; d < polylenmax; d++)
				N_postspikes[d] = 0;
		}

		if (NL_max >= min_group_path)
		{
			int	used[W];
			bool discard = false;
			for (int i = 0; i < W; i++)
			{
				used[i] = 0;
				for (int j = 0; j < N_links; j++)
				{
					if ((links[j, 0] == group[i]) && (links[j, 1] < _numExcitatoryNeurons))
						used[i]++;
				}

				if (used[i] == 1)
					discard = true;
			}

			if (discard == false)
			{
				N_polychronous++;

				// collect together PNG data
				PNGData^ data = gcnew PNGData(N_links, group, t_fired, links);

				// MG changed reporting method from 'write to console' to an event
				FoundPNG(tag, N_polychronous, nnum, N_fired, NL_max, t_fired[N_fired-1], data);
			}
		}
	}

	// Additional comments not provided in the original code have been annotated with my initials (MG)
	// The accuracy of these additional comments has not been established
	void IzhikevichNetwork::SearchCurrentTriplet(String^ tag, int nnum, int npre[W], int &N_polychronous, int polylenmax,
		array<int>^ I_my_pre, array<int>^ D_my_pre,
		array<int>^ group, array<int>^ t_fired, array<int>^ layer, array<int, 2>^ links,
		array<int>^ N_postspikes, array<int, 2>^ I_postspikes, array<int, 2>^ J_postspikes,
		array<int, 2>^ D_postspikes, array<int, 2>^ L_postspikes, array<double, 2>^ C_postspikes, array<double>^ I)
	{
		// This function has been modified from the work of Eugene Izhikevich (see http://www.izhikevich.org/publications/spnet.htm)
		// see BUGFIX for details

		const int latency = _maxDelay;			// maximum latency 

		// MG find the maximum delay within the current set of triplets
		int Dmax = 0;
		for (int i = 0; i < W; i++)
			if (Dmax < D_my_pre[npre[i]])
				Dmax = D_my_pre[npre[i]];

		// MG for each triplet member...
		for (int i = 0; i < W; i++)
		{
			// MG accumulate the neuron index and relative firing time in group[] and t_fired[] arrays
			// MG these can then used to read out the PNG group in event format i.e. firing event == (time, neuron index)
			// MG the layer for each neuron is one more than the layer of its presynaptic neuron in the group (starting from layer 1 for the initial triplet)
			group[i] = I_my_pre[npre[i]];
			t_fired[i] = Dmax - D_my_pre[npre[i]];
			layer[i] = 1;
			
			// MG find common postsynaptic targets (candidates)
			for (int dd = 0; dd < _maxDelay; dd++)									// MG for each delay...
			{
				for (int j = 0; j < _delays_length[group[i], dd]; j++)				// MG and for each offset within the delay band
				{
					int offset = _delays[group[i], dd, j];
					if ((_s[group[i], offset] > C_rel) & (dd >= D_my_pre[npre[i]]))
					{
						// MG found a candidate that has a sufficiently weighty synapse and spikes after i
						// MG store attributes of the candidate indexed by the 'timing' of the postsynaptic spike
						int timing = t_fired[i] + dd + 1;
						J_postspikes[timing, N_postspikes[timing]] = group[i];					// presynaptic
						D_postspikes[timing, N_postspikes[timing]] = dd;						// delay
						C_postspikes[timing, N_postspikes[timing]] = _s[group[i], offset];		// synaptic weight

						// MG the postsynaptic neuron index with the corresponding delay (candidate)
						int p = _post[group[i], offset];
						I_postspikes[timing, N_postspikes[timing]++] = p;						// index of post target
						// MG N_postspikes[timing]: tracks the number of candidates with 'timing'
					}
				}
			}
		}

		// MG initialise membrane potential (v) and recovery vectors (u), and the input vector (I)
		for (int i = 0; i < _totalNeurons; i++)
		{
			_v[i] = -70;
			_u[i] = 0.2 * _v[i];
			I[i] = 0;
		}

		// MG fire all neurons in the current triplet set and their downstream targets
		int N_links = 0;
		int N_fired = W;
		int t_last = _maxDelay + _maxDelay + latency + 1;		// MG latency == D for this algorithm ('latency from the moment spikes arrive to the moment the neuron fires')
		int t = -1;

		while ((++t < t_last) & (N_fired < polylenmax))		// t_last is a variable time limit
		{
			// MG for each candidate, set the input to the synaptic weight
			for (int p = 0; p < N_postspikes[t]; p++) 
				I[I_postspikes[t, p]] += C_postspikes[t, p]; 
 
			// MG update v and u for all neurons and reset the input
			for (int i = 0; i < _totalNeurons; i++)
			{
				_v[i] += 0.5 * ((0.04 * _v[i] + 5) * _v[i] + 140 - _u[i] + I[i]);
				_v[i] += 0.5 * ((0.04 * _v[i] + 5) * _v[i] + 140 - _u[i] + I[i]);
				_u[i] += _a[i] * (0.2 * _v[i] - _u[i]);
				I[i] = 0;
			}

			// MG find neurons that have fired and add them to the candidate list
			for (int i = 0; i < _totalNeurons; i++)
			{
				if (_v[i] >= 30)
				{
					_v[i] = -65;
					_u[i] += _d[i];

					if (N_fired < polylenmax)
					{
						bool enableWarning = true;

						// MG update the firing event list in group[] and t_fired[]
						t_fired[N_fired] = t;
						group[N_fired++] = i;		// MG N_fired tracks the total number of firings

						// MG for each delay and for each offset in the delay band...
						int timing;
						for (int dd = 0; dd < _maxDelay; dd++)
						{
							for (int j = 0; j < _delays_length[i, dd]; j++)
							{
								int offset = _delays[i, dd, j];

								//if ((_s[i, offset] > C_rel) | (i>=_numExcitatoryNeurons))		// MG originally bit-wise OR (delete me)
								if ((_s[i, offset] > C_rel) || (i >= _numExcitatoryNeurons))		// MG if the synapse is weighty or the neuron is inhibitory
								{
									// MG accumulate the attributes of any candidates of fired neuron i at the relevant 'timing'
									timing = t + dd + 1;

									// BUGFIX: it is possible for 'timing' to exceed the maximum index of the _postspikes arrays
									// BUGFIX: it is possible for 'N_postspikes[timing]' to exceed the maximum index of the _postspikes arrays
									// the next lines check these limits
									if ((timing < polylenmax) && (N_postspikes[timing] < _totalNeurons))	// check the limits
									{
										J_postspikes[timing, N_postspikes[timing]] = i;						// presynaptic
										D_postspikes[timing, N_postspikes[timing]] = dd;					// delay
										C_postspikes[timing, N_postspikes[timing]] = _s[i, offset];			// synaptic weight
										I_postspikes[timing, N_postspikes[timing]++] = _post[i, offset];	// index of post target
									}
									else
									{
										if (enableWarning)
										{
											if (timing >= polylenmax)
											{
												RunWarning(String::Format("Warning: timing equalled or exceeded polylenmax: {0}", timing));
												enableWarning = false;
											}

											if (N_postspikes[timing] >= polylenmax)
											{
												RunWarning(String::Format("Warning: N_postspikes[timing] equalled or exceeded polylenmax: {0}", N_postspikes[timing]));
												enableWarning = false;
											}
										}
									}
								}
							}
						}

						// MG push back the time limit (t_last) if necessary
						if (t_last < timing + 1) 
						{
							t_last = timing + 1;
							if (t_last > polylenmax - _maxDelay - 1)
								t_last = polylenmax - _maxDelay - 1;	// MG maximum possible time limit
						}
					}
				}
			}
		}
		
		// MG if the number of fired neurons is over the threshold, calculate links and layers
		if (N_fired > 2 * W)
		{
			N_links = 0;
			int L_max = 0;

			// MG for each fired neuron (other than the original triplet)
			for (int i = W; i < N_fired; i++)
			{
				layer[i] = 0;

				// MG find convergent virtual events that may have been causative in firing neuron i
				for (int p = t_fired[i]; (p > t_fired[i] - latency) && (p >= 0); p--)
				{
					for (int j = 0; j < N_postspikes[p]; j++)
					{
						if ((I_postspikes[p, j] == group[i]) && (J_postspikes[p, j] < _numExcitatoryNeurons)) 
						{
							// MG determine which layer i belongs to (the first layer is 1)
							// MG the layer of i is one greater than the layer of the presynaptic neuron of i
							for (int k = 0; k < i; k++)
							{
								if ((group[k] == J_postspikes[p, j]) && (layer[k] + 1 > layer[i]))
									layer[i] = layer[k] + 1;
							}

							// BUGFIX: it is possible for 'N_links' to exceed the maximum index of the links arrays
							// the next line checks this limit
							if (N_links < 2 * W * polylenmax)	// check the limit
							{
								// MG add attributes to the links[,] array
								links[N_links, 0] = J_postspikes[p, j];		// MG pre
								links[N_links, 1] = I_postspikes[p, j];		// MG post
								links[N_links, 2] = D_postspikes[p, j];		// MG delay
								links[N_links++, 3] = layer[i];				// MG layer
								if (L_max < layer[i])
									L_max = layer[i];		// MG track the maximum layer in L_max
							}
						}
					}
				}
			}
			
			// MG discard the group if any member of the original triplet has no links to excitatory neurons
			bool discard = false;
			int used[W];
			for (int i = 0; i < W; i++)
			{
				used[i] = 0;
				for (int j = 0; j < N_links; j++)
					if ((links[j, 0] == group[i]) && (links[j, 1] < _numExcitatoryNeurons))
						used[i]++;
				if (used[i] == 1)
					discard = true;
			}

			// MG if the maximum layer reaches a threshold...
			if ((discard == false) && (L_max >= min_group_path))
			{
				// MG we have found a polychronous group
				// MG accumulate the count and report
				N_polychronous++;

				// collect together PNG data
				PNGData^ data = gcnew PNGData(N_links, group, t_fired, links);

				// MG changed reporting method from 'write to console' to an event
				FoundPNG(tag, N_polychronous, nnum, N_fired, L_max, t_fired[N_fired-1], data);
			}
		}

		// MG reinitialise the values of N_postspikes in the utilised range to zero...
  		for (int dd = Dmax; dd < t_last; dd++)
			N_postspikes[dd] = 0;

		// MG ...and any additional values where the time limit has been pushed back
		if (t_last == polylenmax - _maxDelay)
		{
			for (int dd = t_last; dd < polylenmax; dd++)
				N_postspikes[dd] = 0;
		}
	}
}